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Abstract 

We describe a simple coarse-grained model which is suited to study lipid layers and their phase transitions. Lipids 
are modeled by short semiflexible chains of beads with a solvophilic head and a solvophobic tail component. They are 
forced to self- assemble into bilayers by a computationally cheap 'phantom solvent' environment. The model reproduces 
the most important phases and phase transitions of monolayers and bilayers. Technical issues such as Monte Carlo 
parallelization schemes are briefly discussed. 
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1. Introduction 

Lipid bilayers are the main components of biolog- 
ical membranes and omnipresent in all living mat- 
ter (0). At high temperatures bilayers assume the 
so-called 'liquid' state (I/a), where lipids are highly 
mobile and have many chain defects. In nature, this 
state is the most frequent. If one decreases the tem- 
perature, pure one-component lipid bilayers undergo 
a prominent phase transition, the 'main' transition, 
which is characterized by dropping lipid mobility, 
dropping number of chain defects, and dropping area 
per lipid. The structure of the low temperature 'gel' 
phase depends on the bulkiness and interactions of 
the head groups. For small head groups, the chains 
are oriented normal to the bilayer (L/3 phase), for 
larger head groups, they are tilted {Lp'). In the latter 
case, the main transition occurs in two steps, and an 
undulated intermediate phase emerges, the 'ripple' 
phase P/3/ . If head groups are large and weakly inter- 
acting, such as ether-linked phosphatidylcholines, 
the system assumes a phase where both oppos- 
ing lipid layers are fully interdigitated (2). 



In this paper, we present a lipid model which is 
suited for studying lipid bilayers. We will first apply 
it to lipid monolayers (Sec. 2) and show that it re- 
produces the generic features of fatty acid monolay- 
ers. Then we introduce an environment model which 
forces the model lipids to self-assemble into bilay- 
ers, and discuss the resulting bilayer phases (Sec. 3). 
Selected technical issues regarding the Monte Carlo 
implementation are discussed in the Appendix. 

2. Lipids and Monolayers 

The lipids are represented by chains of n — 1 'tail' 
beads with diameter cr^, attached to one 'head' bead 
with diameter ah- Beads that are not direct neigh- 
bors along their chain interact with a truncated and 
shifted Lennard- Jones potential. 
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VLj{r) = otherwise, with Vc chosen such that 
VLj{r) is continuous at r = i?c- The parameter a is 
the arithmetic mean of the diameters of the two in- 
teracting beads. Head-head interactions and head- 
tail interactions are purely repulsive, which is en- 
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sured by choosing Rq = a. Tail-tail interactions have 
an attractive part, Rq = 2(j. 

Within a chain, beads are connected by bonds 
of length d subject to the weakly nonlinear spring 
potential (FENE potential) 
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for \d-do\ < ds (2) 



and Vs{d) = oo otherwise, where do is the equi- 
librium spring length, kg the spring constant, and 
the logarithmic cutoff ensures that the spring length 
never exceeds do -\- ds. 

In addition, a bending potential 



Va=ka{l-COsO) 



(3) 



is imposed, which acts on the angle between sub- 
sequent bonds. 

The parameters at and e provide 'natural' length 
and energy units of the system. In these units, we 
use ks = lOOe (very stiff bonds), do = 0.7 at, ds = 
0.2crt, ka = 4.7e. The values are motivated by simple 
considerations that map our chains on hydrocarbon 
chains (0;0); the 'matching' should not be taken too 
literally, since the model is not designed to describe 
experiments on a quantitative level. The size of the 
head group, ah, and the chain length, n, are model 
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Fig. 1. Monolayer phase diagrams, (a) Generic phase diagram 
for fatty acid monolayers (after Ref. (3)). LE is the liquid 
expanded phase, the other phases are ordered hexatic liquids. 
The chains are untilted in LS, and they tilt in different 
directions in L2 and Ov. (b) Phase diagram of our lipid 
model. From Ref. (0). 



parameters that allow to study the influence of the 
head group bulkiness and the chain length on the 
phase behavior (I5I). Unless stated otherwise, they 
are chosen ah — 1.1 and n = 7. 

To evaluate the properties of the lipid model, we 
first consider monolayers of lipid at an air- water in- 
terface. Such monolayers have been studied for a 
long time as experimentally accessible model sys- 
tems for lipid layers (0; 0) . The monolayer equiva- 
lent of the main transition is a transition encoun- 
tered upon compression of the monolayer from a 'liq- 
uid expanded' (LE) phase to a 'liquid condensed' 
phase. As in bilayers, the ordered 'liquid condensed' 
phase exists in several modifications, which differ, 
among other things, in the tilt order of the chains 
(L2, LS, or Ov phase, see Fig. la). In the simulations, 
the water surface can be replaced by suitable exter- 
nal potentials. With smooth harmonic potentials of 
width ~ cr, we obtain the phase diagram shown in 
Fig. lb) ( 4). It is in good qualitative agreement with 
the experimental phase diagram. Fig. la). 



3. Phantom solvent and self-assembly 

Having formulated a reasonable lipid model, we 
must now force the 'lipids' to self-assemble into bi- 
layers. In nature, self-assembly is caused by the in- 
teraction with the surrounding water, hence we must 
add an appropriate representation for the environ- 
ment. This is done by introducing a recently pro- 
posed, simple and very efficient environment model: 
The 'phantom solvent' model (6). Explicit 'solvent' 
particles are added to the system, which however do 
not interact with each other, only with lipid beads 
(by means of repulsive interactions, Eq. (1)). Phys- 
ically, the solvent probes the accessible free volume 
in the presence of lipids on the length scale of the 
solvent diameter as. Therefore, it promotes lipid ag- 
gregation, and the lipids self-assemble to bilayers 
(see Fig. 2). Compared to other explicit solvent mod- 
els |3), the phantom solvent environment has the 
advantage of having no internal structure, it thus 
transmits no indirect interactions between different 
bilayer regions and/or periodic images of bilayers. 
Furthermore, it is cheap - in Monte Carlo simula- 
tions, less than 10 % of the computing time is spent 
on the solvent. Compared to implicit solvent mod- 
els (8), where the solvent is replaced by effective 
lipid interactions, it has the advantage that no tun- 
ing of potentials is required, and it can also be used 
to study solvent dynamics. For example, with DPD 
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Fig. 2. Snapshot of lipid bilayers a) fluid bilayer L^, b) tilted 
gel L^/ c) asymmetric ripple P^/ 

dynamics, one can study the effect of hydrodynamic 
coupling between membranes and the surrounding 
fluid. 

In our work, the solvent diameter was chosen ag = 
(Jh- Single head beads are soluble, (i.e., they do not 
demix with solvent) if the free solvent density is 
less than pfree ^ 2.6/<j^. At sufficiently low tem- 
peratures, the lipids self-assemble into bilayers (see 
Fig. 2). The properties of these bilayers will be dis- 
cussed in detail elsewhere Here we just cite 
some of the main results. Like the monolayers, the 
bilayers exhibit a main transition. For small heads 
{ah = 0.9cr), the gel phase is untilted, ie., we ob- 
tain an La phase. For larger heads {ah = l.lcr), the 
structure of the gel phase depends on the free sol- 
vent density Pfree- We note that the solvent entrop- 
ically penalizes lipid/solvent interfaces and thus ef- 
fectively creates an attractive depletion interaction 
between the beads next to these interface, ie., the 
head beads. The strength of this interaction is pro- 
portional to Pfree- At loW /Ofree, thc gcl phaSC IS lU- 

terdigitated {L'^^)^ at moderate pfree > 1.2/cr^, it is 
tilted {Lp'). Hence weak head attraction leads to the 
formation of the interdigitated phase, and moderate 
head attraction to the tilted gel phase, in agreement 
with experiments. 

Most rewardingly, we also recover the ripple 
phase Pf3' which intrudes between the tilted gel 
phase and the ffuid phase. A snapshot is shown in 
Fig. 2 c). The two main experimental rippled states, 
the 'asymmetric' and the 'symmetric' rippled state, 
are recovered in simulations, with properties that 
are very similar to experimental properties (9). A 
similar structure than that of our asymmetric ripple 
has been found recently in a (much more involved) 
atomistic simulation of a Lecithine bilayer (11). Our 
simulations show that this structure is generic, in 
the sense that it can be reproduced with a coarse- 
grained model, and that it is closely related to the 
structure of the corresponding symmetric rippled 




Fig. 3. Fluid membrane with two embedded coarse-grained 
transmembrane proteins 

state (191). 



4. Conclusions and Outlook 

To conclude, we have presented a versatile coarse- 
grained model that allows to study lipid monolayers 
and self-assembled bilayers and reproduces their 
most important internal phase transitions. It can 
be used to study a variety of questions related to 
membrane biophysics where atomic details do not 
matter, but the characteristic molecular features 
of lipids are still important. For example, we are 
currently applying it to study lipid-mediated inter- 
action mechanisms between proteins An example 
snapshot is shown in Fig. 3. 
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Appendix: Technical Remarks 

The Monte Carlo simulations described above 
were carried out at constant pressure P with pe- 
riodic boundary conditions in a simulation box of 
variable size and shape: The simulation box is a 
parallelepiped spanned by the vectors (La^,0,0), 
(5iL^,L^,0), and {s2Lx, ssLy, L^). All and Si 
are allowed to fluctuate. In addition, it is sometimes 
convenient to work in a semi-grand canonical en- 
semble with fluctuating number Ng of solvent beads 
and given solvent chemical potential jUs- Hence 
we have three types of possible trial Monte Carlo 
moves: 

- Moves that change the positions of beads. 

- Moves that changes the volume and/or shape of 
the box: Random increments drawn randomly 
from a symmetric distribution with mean zero 
are added to or Si. All bead coordinates are 
rescaled accordingly. 
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- Moves that changes the number Ng of solvent 
particles: We first decide with probability 1/2 
whether to attempt a solvent removal or a solvent 
addition. Then, we choose randomly the particle 
to be removed, or the position of the particle to 
be added. 

The moves are accepted or rejected according to 
one of the standard Monte Carlo schemes (e.g.. 
Metropolis), with the effective Hamiltonian (Il2l) 

ifeff = H+PV-^,N,-kBTln [{V/Vo)''/N,\], (4) 

where H is the interaction energy, V = L^LyLz 
the volume of the simulation box, Vq an arbitrary 
reference volume (e.^., Vo = a^), and the total 
number of beads (solvent and lipids). The La are not 
allowed to fall below a given threshold, otherwise 
the move is rejected. 

We close with a remark on parallelization. For 
large scale applications, our Monte Carlo code has 
been parallelized geometrically. One commonly used 
spatial decomposition scheme for systems with short 
range interactions (see, e.g., the review (0)) pro- 
ceeds as follows: The simulation box is divided into 
domains which are distributed on the processors. 
These are further subdivided into labelled subdo- 
mains such that subdomains with the same label are 
separated by a distance larger than the maximum 
interaction range. Subdomains with the same label 
are then processed in parallel. This algorithm is rel- 
atively straightforward, yet it has the drawback that 
it does not strictly fulfill detailed balance: Within 
a move for given subdomain label a, particles can 
cross a subdomain boundary only in one direction 
{i.e., leaving the set of subdomains a). If the differ- 
ent sets a are processed equally often, the final dis- 
tribution is presumably not affected. Nevertheless, 
we feel uncomfortable with this method and favor 
a variant of a parallelization scheme recently pro- 
posed by Uhlherr et al. (14). The idea is to define 
'active regions' and assign them to different proces- 
sors. A possible decomposition scheme is shown in 
Fig. 4. Only particles with centers in the active re- 
gions are moved, and moves that take a particle out- 
side of its active region are rejected. To ensure that 
the algorithm remains ergodic, the active regions are 
periodically redefined. One easily checks that indi- 
vidual bead moves fulfill detailed balance. We note 
that the active regions must not necessarily have the 
same size, and furthermore moves in interesting re- 
gions can be more frequent. This feature makes the 
concept of 'active regions' interesting even for appli- 
cations on scalar computers. 





Fig. 4. Sketch of the domain decomposition scheme used for 
the parallehzation of the Monte Carlo algorithm (a variant 
of Uhlherr et al. 13) V The system is covered with a grid 
of active regions (shaded rectangles, filled particles). The 
distance between active regions must exceed the interaction 
range between beads. The offset of the active grid (arrow) 
changes periodically and is chosen randomly. 
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